﻿/* Copyright 2007-2009 dnAnalytics Project.
 *
 * Contributors to this file:
 * Jurgen Van Gael
 * 
 * Redistribution and use in source and binary forms, with or without modification,
 * are permitted provided that the following conditions are met:
 * 
 * * Redistributions of source code must retain the above copyright notice, this 
 *   list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above copyright notice, 
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 * * Neither the name of the dnAnalytics Project nor the names of its contributors
 *   may be used to endorse or promote products derived from this software without
 *   specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
 * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
using System;

namespace dnAnalytics.Mcmc
{
    /// <summary>
    /// Slice sampling produces samples from distribition P by uniformly sampling from under the pdf of P using
    /// a technique described in "Slice Sampling", R. Neal, 2003. All densities are required to be in log space.
    /// 
    /// The slice sampler is a stateful sampler. It keeps track of where it currently is in the domain
    /// of the distribution P.
    /// </summary>
    public class UnivariateSliceSampler : McmcSampler<double>
    {
        /// <summary>
        /// Evaluates the log density function of the target distribution.
        /// </summary>
        private readonly DensityLn<double> mPdfLnP;
        /// <summary>
        /// The current location of the sampler.
        /// </summary>
        private double mCurrent;
        /// <summary>
        /// The log density at the current location.
        /// </summary>
        private double mCurrentDensityLn;
        /// <summary>
        /// The number of burn iterations between two samples.
        /// </summary>
        private int mBurnInterval;
        /// <summary>
        /// The scale of the slice sampler.
        /// </summary>
        private double mScale;

        /// <summary>
        /// Constructs a new Slice sampler using the default <see cref="System.Random"/> random 
        /// number generator. The burn interval will be set to 0.
        /// </summary>
        /// <param name="x0">The initial sample.</param>
        /// <param name="pdfLnP">The density of the distribution we want to sample from.</param>
        /// <param name="scale">The scale factor of the slice sampler.</param>
        /// <exception cref="ArgumentOutOfRangeException">When the scale of the slice sampler is not positive.</exception>
        public UnivariateSliceSampler(double x0, DensityLn<double> pdfLnP, double scale) :
            this(x0, pdfLnP, 0, scale)
        {
        }

        /// <summary>
        /// Constructs a new slice sampler using the default <see cref="System.Random"/> random number generator. It 
        /// will set the number of burnInterval iterations and run a burnInterval phase.
        /// </summary>
        /// <param name="x0">The initial sample.</param>
        /// <param name="pdfLnP">The density of the distribution we want to sample from.</param>
        /// <param name="burnInterval">The number of iterations in between returning samples.</param>
        /// <param name="scale">The scale factor of the slice sampler.</param>
        /// <exception cref="ArgumentOutOfRangeException">When the number of burnInterval iteration is negative.</exception>
        /// <exception cref="ArgumentOutOfRangeException">When the scale of the slice sampler is not positive.</exception>
        public UnivariateSliceSampler(double x0, DensityLn<double> pdfLnP, int burnInterval, double scale)
        {
            mCurrent = x0;
            mCurrentDensityLn = pdfLnP(x0);
            mPdfLnP = pdfLnP;
            Scale = scale;
            BurnInterval = burnInterval;

            Burn(BurnInterval);
        }

        /// <summary>
        /// Gets or sets the number of iterations in between returning samples.
        /// </summary>
        /// <exception cref="ArgumentOutOfRangeException">When burn interval is negative.</exception>
        public int BurnInterval
        {
            get { return mBurnInterval; }

            set
            {
                if (value < 0)
                {
                    throw new ArgumentOutOfRangeException(Properties.Resources.ParameterCannotBeNegative);
                }
                    mBurnInterval = value;
            }
        }

        /// <summary>
        /// Gets or sets the scale of the slice sampler.
        /// </summary>
        public double Scale
        {
            get { return mScale; }

            set
            {
                if (value <= 0.0)
                {
                    throw new ArgumentOutOfRangeException(Properties.Resources.NotPositive);
                }
                    mScale = value;
            }
        }

        /// <summary>
        /// This method runs the sampler for a number of iterations without returning a sample
        /// </summary>
        private void Burn(int n)
        {
            for (int i = 0; i < n; i++)
            {
                double x_l = mCurrent;
                double x_r = mCurrent;
                double xnew = mCurrent;

                // The logarithm of the slice height.
                double lu = System.Math.Log(RandomNumberGenerator.NextDouble()) + mCurrentDensityLn;
                
                // Create a horizontal interval (x_l, x_r) enclosing x.
                double r = RandomNumberGenerator.NextDouble();
                x_l = mCurrent - r * Scale;
                x_r = mCurrent + (1.0 - r) * Scale;
                
                // Stepping out procedure.
                while (mPdfLnP(x_l) > lu) { x_l -= Scale; }
                while (mPdfLnP(x_r) > lu) { x_r += Scale; }

                // Shrinking: propose new x and shrink interval until good one found.
                while (true)
                {
                    xnew = RandomNumberGenerator.NextDouble() * (x_r - x_l) + x_l;
                    mCurrentDensityLn = mPdfLnP(xnew);
                    if (mCurrentDensityLn > lu)
                    {
                        mCurrent = xnew;
                        mAccepts++;
                        mSamples++;
                        break;
                    }
                    if (xnew > mCurrent)
                    {
                        x_r = xnew;
                    }
                    else
                    {
                        x_l = xnew;
                    }
                }
            }
        }

        /// <summary>
        /// Returns a sample from the distribution P.
        /// </summary>
        public override double Sample()
        {
            Burn(BurnInterval + 1);

            return mCurrent;
        }
    }
}
